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ABSTRACT 

We present the results of fully 3-D hydrodynamic simulations of the gravitational collapse of isolated, 
turbulent molecular cloud cores. Starting from initial states of hydrostatic equilibrium, we follow the 
collapse of both singular and nonsingular logatropic cores until the central protostar has accreted > 90% 
of the total available mass. We find that, in the collapse of a singular core with access to a finite mass 
reservoir, the mass of the central protostar increases as M^cc oc until it has accreted ^ 35% of the 
total available mass. For nonsingular cores of fiducial masses 1, 2.5, and 5 Mq, we find that protostellar 
accretion proceeds slowly prior to the formation of a singular density profile. Immediately thereafter, 
the accretion rate in each case increases to ~ 10~^ Mq yr~^, for cores with central temperature Tc = 10 
K and truncation pressure Pg = 1.3 x lO^fcs cm~^ K. It remains at that level until half the available 
mass has been accreted. After this point, the accretion rate falls steadily as the remaining material is 
accreted onto the growing protostellar core. We suggest that this general behaviour of the protostellar 
accretion rate may be indicative of evolution from the Class to the Class I protostellar phase. 

Subject headings: stars: formation — methods: numerical — equation of state 



1. INTRODUCTION 

A crucial prerequisite for any comprehensive theory of 
isolated star formation is an understanding of the gravi- 
tational collapse of molecular cloud cores^. Larson (1969) 
and Penston (1969) pioneered the study of the collapse of 
the isothermal sphere from both analytical and numerical 
perspectives. Subsequently, Shu (1977) defined the long- 
standing paradigm for isolated low mass star formation — 
the "inside-out" or "expansion wave" collapse — with his 
elegant self-similar solution for the collapse of a singular 
isothermal sphere (SIS). 

A common feature of all three of the aforementioned 
studies, and most star-formation studies which have fol- 
lowed, is their adoption of the isothermal equation of state 
(EOS). Recent evidence supports the assertion that some 
observed molecular cores are well-fit by isothermal mod- 
els. For example, Alves, Lada, & Lada (2001) demon- 
strated convincingly that the density distribution of the 
dark cloud Barnard 68 is very well fit by a critically stable 
Bonnor-Ebert sphere (Bonnor 1956; Ebert 1955). 

There is, however, mounting evidence which indicates 
that more massive molecular cores possess properties 
which cannot be explained by the isothermal EOS. Sur- 
veys of prestellar cores in a variety of star-forming com- 
plexes show that massive cores have strong nonthermal 
linewidths which are greater than their thermal linewidths 
(Caselli and Myers 1995). Furthermore, it has been shown 
by Fuller & Myers (1992) that even low mass cores have 
significant nonthermal linewidths. 

The SIS model for molecular cloud core collapse pre- 



dicts a time-invariant rate of accretion onto the central 
protostar (Shu 1977), which gives rise to a large asso- 
ciated spread between the formation times of high and 
low mass stars. In the high-mass star formation regime, 
this is problematic: the constant accretion rate predicted 
by SIS models is too low to form massive stars on the 
timescales suggested by observations (Caselli and Myers 
1995; McLaughlin & Pudritz 1997, MP97 hereafter). Fos- 
ter & Chevaher (1993) and Boss and Black (1982) have 
shown that it is possible to obtain variable accretion rates 
from isothermal collapse models if, for example, the mass 
reservoir is finite or there are initial deviations from the 
SIS density profile. It is not clear, however, that the mass 
accretion rates generated by these models vary in a way 
which is consistent with observations. The likely conclu- 
sion is that isothermal models only apply to the formation 
of low mass stars. 

The SIS collapse model assumes hydrostatic initial con- 
ditions; it has been suggested that such a singular con- 
figuration could arise through the subsonic collapse of a 
Bonnor-Ebert sphere (Shu, Adams, & Lizano 1987). Sim- 
ulations of the collapse of critical Bonnor-Ebert spheres, 
by contrast, indicate that, as the density achieves a singu- 
lar form (i.e. when p{r) takes on an profile), 44% of the 
mass is infalling at a few times the sound speed (Foster & 
Chevalier 1993; FC93 hereafter). Such high infall velocities 
have not yet been detected in collapsing cores, but whether 
this is because they do not exist, or because their detection 
is beyond the limits of current observational techniques is 
not yet known (see FC93 for discussion). 

One proposed alternative to the isothermal EOS is the 



^ To minimize ambiguity, we will use the terms "cloud" , "clump" , and "core" in the following hierarchical manner; molecular clouds contain 
dense condensations called clumps, in which entire star clusters may form; clumps, in turn, contain cores, from which single or binary stars 
may form. The centre of a core shall be denoted herein with the terms "central object" or "central protostar". 
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logatropic equation of state, so named for its logarith- 
mic dependence of pressure on density. First proposed 
by Lizano and Shu (1989), the logatrope is an empirically- 
motivated EOS which attempts to account for the pres- 
ence of nonthermal or turbulent pressure within molecular 
cores. In the so-called "mixed" form of Lizano and Shu 
(1989), the logatrope was essentially just the isothermal 
equation of state with a correction added to account for 
turbulent support. McLaughHn & Pudritz (1996) (MP96 
hereafter) proposed the "pure" form of the logatrope, 
which eliminated the linear (isothermal) dependence of 
pressure on density, leaving only the logarithmic depen- 
dence. In this form, the logatrope can describe the prop- 
erties of both low and high mass cores (we refer here pri- 
marily to the properties of molecular cores as determined 
by spatially unresolved measurements, such as single mea- 
surements of the thermal and nonthermal components of 
a core's velocity dispersion). While nonthermal support is 
certainly important in high mass cores, low mass cores are 
also characterized by some degree of nonthermal support. 

The logatrope is capable of addressing some of the dif- 
ficulties with the isothermal models. For example, it 
is specifically formulated to account for the nonthermal 
linewidths of observed cores. Also, because they predict 
a variable protostellar mass accretion rate {Mace oc t^), 
self-similar models for the collapse of a singular logatropic 
sphere (SLS) give rise to a relatively small spread between 
the formation times of high and low mass stars. 

The general consensus at this time favours models which 
produce subsonic molecular cloud core collapse. Because 
the logatrope is a softer EOS than the isothermal EOS (i.e. 
P depends less than linearly on p), it has been suggested 
that it might allow for a gentler, subsonic approach to 
the singular density profile (MP96). This contrasts with 
the strongly supersonic flow found in simulations of the 
collapse of Bonnor-Ebert spheres (FC93). The logatrope 
has the potential to produce gentler collapses because it 
accounts for some measure of turbulent and magnetic sup- 
port within a molecular cloud core. 

There is direct observational evidence that the internal 
(i.e. spatially resolved) structure of individual turbulent 
cores matches the predictions of logatropic models. For 
example, van der Tak et al. (2000) examined the structure 
of the envelopes around 14 massive young stars and found 
that the density profiles of these envelopes have power law 
distributions, n cx r", with indices a = 1.0 — 1.5. This 
is consistent with the logatropic model, which predicts 
a = 1.0 for cores in equilibrium and a = 1.5 for collaps- 
ing cores. Further, these values of a indicate density pro- 
files which are significantly fiatter than the a = 2.0 ± 0.3 
commonly found for low mass cores, and are inconsistent 
with the value of a = 2.0 predicted by the SIS model 
for an equilibrium isothermal core. These results indicate 
that nonthermal support mechanisms dominate in mas- 
sive yoking stellar objects, while thermal pressure dom- 
inates in objects of lower mass. Similarly, Henning ct 
al. (1998) and Colome, di Francesco, & Harvey (1996) 
find a = 0.75—1.5 for the envelopes around massive 
Herbig Ae/Be stars. Other evidence is provided by Os- 
orio, Lizano, & D'Alessio (1999), who examined the dust 
thermal spectra of several hot molecular cores, which are 
thought to be the sites of massive star formation. They 



showed that the thermal spectra of these cores were best 
fit by models whose envelopes had the density profile of 
collapsing logatropic spheres. 

The equilibria, stability criteria, and self-similar collapse 
solutions for pure logatropic spheres have been derived by 
MP96 and MP97. They showed that logatropic spheres 
have two equilibria: one is a stable equilibrium with finite 
central density and a p oc envelope, and the other is 
an unstable equilibrium with a singular (p oc r~^) density 
profile throughout. These crjuilibria arc roughly analogoiis 
to those of the isothermal sphere, namely the hydrostatic 
Bonnor-Ebert sphere and the unstable SIS. MP97 also de- 
rived self-similar collapse solutions for the SLS, similar to 
those developed by Shu (1977) for the SIS. 

In this paper, we extend the work of MP97 by perform- 
ing fully three-dimensional numerical hydrodynamic sim- 
ulations of the gravitational collapse of logatropic spheres, 
both singular and nonsingular. Our goal is to develop a set 
of numerical collapse solutions for the logatrope which is 
analogous to those developed for the Bonnor-Ebert sphere 
by FC93. We present results of simulations which de- 
scribe the approach to a singular density profile in the 
collapse of an initially hydrostatic, nonsingular logatropic 
sphere. Because the collapse of a nonsingular logatrope is 
not suited to exploration by the usual self-similar analyt- 
ical techniques, we are motivated to use a more flexible 
numerical approach. 

In §2, we review the analytic theory of logatropic; equi- 
libria and collapses. Section 3 describes our numerical 
method and discusses our testing procedures. Sections 4 
and 5 present the results of our simulations of the singular 
and nonsingular logatropic collapses (scaled such that the 
total core mass is 1 Mq), respectively. Section 6 discusses 
the results of simulations of the collapse of higher-mass 
nonsingular cores and relates them to the observed prop- 
erties of protostellar cores. We conclude with discussion 
and a summary in sections 7 and 8. 

2. LOGATROPIC COLLAPSE THEORY 

2.1. The Logatropic Equation of State 

The logatrope was proposed by Lizano and Shu (1989) 
as a way to model the observed nonthermal component 
of the velocity dispersion in molecular clouds. Their ap- 
proach was to separate the equation of state into thermal 
and nonthermal components, where they assumed that 
the nonthermal component was turbulent in origin, hence 
-Ptot = -Puicrm + -Piurb- In this formulation, -Puicrm is iden- 
tified with the isothermal contribution to the pressure, 
-Ptherm = ^iso = c^p. The form of Pturb is chosen to fit 
the observed correlation between the nonthermal velocity 
dispersion and density in molecular clouds; 
HP 

^ = (A^;nt)' cx (1) 
dp 

Pturb (X In ( ) , (2) 

VPref/ 

where Pref is the density at some reference point. Equa- 
tion (1) follows simply from the empirical scalings re- 
ported by Larson (1981), namely that the nonthermal ve- 
locity dispersion and density within molecular cores scale 
as Awnt oc r^/^ and p oc r~^. MP96 noted that this 
isothermal/logarithmic "mixed" equation of state had the 
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property that P/ p decreased with radius, contrary to ob- 
servations. In its place, MP96 proposed a "pure" form of 
the logatrope: 



= 1 + A In 



(3) 



where the subscript "c" denotes "central" values and Pc = 
(tIPc- Here (T^ = kBTc/ pmn is the central velocity dis- 
persion and Tc is the central temperature of the molecular 
core. To account for the fact that the linewidths in the cen- 
tres of molecular clouds are observed to be purely thermal, 
equation (3) is formulated such that, when p = p^, it gives 
P = Pc = (J^Pc- Note that it is only at the precise centre of 
the molecular cloud core that this correspondence between 
the isothermal and logatropic models is established — there 
is no extended central isothermal region. 

The constant A in equation (3) is an adjustable param- 
eter, and it is always positive. The best-fit value of A 
was determined by MP96, who fit equation (3) to observa- 
tional data on both low mass (Fuller & Myers 1992) and 
high mass (Caselli and Myers 1995) GMC cores, obtaining 
A = 0.20 ± 0.02. In this work, we take A = 0.2. 

It is the pure logatropic equation of state of MP96 whose 
properties we explore in this paper. Henceforth, equation 
(3) shall be referred to simply as "the logatrope" , with the 
understanding that no further discussion shall be made of 
other logarithmic equations of state. 

2.2. Logatropic Equilibria 

In order to examine the equilibria of logatropic spheres, 
we begin with the equations governing spherically- 
symmetric flow in a self-gravitating fluid: 



du du dcp 1 dP Q 
dt dr dr p dr 



dp 
dt 



--A'r'f^\=^^Gp 



(4) 



(5) 



(6) 



where u{r, t) is the fluid velocity at radius r and time t, 
p(r,t) is the density, P{r,t) is the pressure, and 4){r,t) is 
the gravitational potential. 

The equilibria of the logatrope can be found by setting 
terms in u and d/dt to zero in equations (4)-(6), insert- 
ing the logatropic EOS, and solving. One solution that 
emerges has a singular density proflle, of the form 



Pc M 9 \r) ' 



(7) 



where ro, the scale radius of the core, is deflned as Tq = 
Qa'l/A'KGpc- This is the logatropic analogue of the sin- 
gular isothermal sphere. The second solution, which can 
be found by numerical integration, represents a pressure- 
truncated sphere with a finite central density; this is the 
logatropic analogue of the Bonnor-Ebert sphere. The de- 
tailed properties of the two solutions are developed in 
MP96, but their density profiles are reproduced here for 
the reader's convenience (see Fig. 1). We note that the 
nonsingular profile follows the singular profile closely over 
a broad range of radii. 



MP96 determined the stability criteria for both the sin- 
gular and nonsingular equilibria. They found that singu- 
lar logatropcs of any size are unstable to collapse, but that 
nonsingular logatropcs of radius R are stable if 



R < i?crit — fO 




exp 



(8) 



9 4, 

The critical radius of an ^ = 0.2 pressure-truncated loga- 
tropic sphere is i?crit/''o = 24.37. 



2.3. Logatropic Collapse 

In §4, we compare our SLS simulation results to the 
analytic collapse solutions of MP97. We refer the inter- 
ested reader to that paper for a complete treatment of the 
matter, but we highlight some of the relevant results here. 

The SLS collapse solution is qualitatively similar to 
the SIS collapse: both are characterized by an expan- 
sion wave which begins at the origin and moves outward 
through a static medium, causing the material it passes 
to collapse inward onto a forming protostar. In the case 
of the SLS, the radius of the expansion wave grows as 
rew = {APc'KG/Sy^'^t'^. The passage of the expansion 
wave causes material at radii r < rew to collapse, while 
material at r > rew remains stationary. The expansion 
wave reaches the edge of the core at time 

tsuri = 4V2/7rfff ~ l.Sis , (9) 

where tg = (37r/32pG)^/^ is the mean free-fall time in the 
core. 

The asymptotic collapse solutions for small radii or late 
times are the same for both the SIS and SLS collapse mod- 
els: p oc r~'^/^ and v oc — r~^/-^, where v is the mean veloc- 
ity of the bulk flow. One crucial diff'erence between the two 
models is their widely-differing mass accretion rates. The 
SIS collapse is characterized by a constant rate of accre- 
tion onto the central object, M = 0.975a^/G, where a is 
the isothermal sound speed. In contrast, the SLS collapse 
exhibits a time-varying accretion rate, 

M(0,t) = iiAPc^nGf/^"^ . (10) 

Here, mo is a dimensionless parameter — the reduced mass 
of the collapsed central object — arising from the self- 
similar collapse analysis. The primary function of TOq is to 
determine the density normalization of a given SLS, and 
thereby to set the mass accretion rate. For singular loga- 
tropcs. Too = 6.67 X 10~^, which is the value used herein. 
This scaling of the mass of the protostellar object with 
time provides one convenient test of the accuracy of our 
simulations. 

3. NUMERICAL METHODS 

3.1. Motivation for Numerical Approach 

We adopt the numerical approach as a means of ex- 
tending in time the self-similar SLS collapse solutions of 
MP97, and also as a tool for exploring the collapse of ob- 
jects which do not lend themselves to the self-similar ap- 
proach, such as the nonsingular logatrope. Further, the 
development of a numerical technique is preparatory to 
future work wherein we will generalize our approach to 
study the collapse of rotating, magnetized logatropic cores. 
This purely hydrodynamic study serves to test our ability 




Fig. 1. — Density profiles for equilibrium logatropic spheres with A=0.20. The nonsingular and singular solutions are represented by the 
solid and dashed lines, respectively. The vertical line represents the truncation radius of a critically stable cloud (with R/ro = 24.37). 



to perform accurate 3-D logatropic collapses at feasible 
resolutions. Magnetized, rotating collapses are likely to 
be characterized by significant deviations from spherical 

symmetry, so their study will be greatly facilitated by the 
use of 3-D magnetohydrodynamic (MHD) simulations. 

3.2. The ZEUS-MP Code 

Our simulations were conducted using ZEUS-MP (Stone 

& Norman 1992a, b), a parallel MHD code with a built-in 
gravity solver, which we obtained from the National Center 
for Supercomputing Applications (NCSA) at the Univer- 
sity of Illinois. ZEUS-MP uses a finite-difference Eulerian 
grid scheme to solve the hydrodynamic equations and a 
multigrid algorithm to solve the Poisson equation for self- 
gravity. The algorithms are formally second order in both 
time and space. 

The simulations were run on the McMastcr component 
of SHARC-NET, a distributed Beowulf cluster shared by 
several Ontario universities. McMaster's component is an 
AlphaServer SC with 112 processors and 1 GB of RAM 
per processor. 

Simulations wore conducted on uniform Cartesian grids, 
with resolutions of both 127'^ and 255^ cells. The use of a 
uniform Cartesian grid was dictated in part by the ZEUS- 
MP gravity solver, which could not be made to perform 
accurately on other grid types. 

3.3. Modifications to the Code 

We made several modifications to the NCSA version of 
ZEUS-MP. The calculation of the pressure term in equa- 
tion (4) was modified to accommodate the logatropic EOS. 
ZEUS-MP updates the velocities in each coordinate direc- 
tion separately (Stone & Norman 1992a). Thus, in the 
a;-direction, the change in velocity due to the pressure gra- 
dient is typically computed as: 

v^,i{t + M) - v^,i{t) _ P^{t) - P,_i(0 

At Aa;,(p,(t)+p,_i(i))/2 ^ ' 

where v^^iit) is the velocity in the a;-direction at grid 
point i and time t, Pi{t) is the pressure, piit) is the den- 
sity. At is the time step, and Axj is the grid spacing in 



the a;-direction. ZEUS-MP uses a staggered grid method, 
whereby the grid of vector quantities is offset from that 
of scalar quantities by half a grid spacing. Thus, the ve- 
locities in equation (11) are computed at a point offset by 
Aa;j/2 from the densities and pressures. Equation (11) as- 
sumes that the pressures have been calculated previously 
according to one of the equations of state built into ZEUS- 
MP. For a given equation of state, however, we can calcu- 
late the change in velocity directly, without first calculat- 
ing the pressures separately. Casting the logatropic EOS 
into dimensionless units {P = 1 + A\np), we can write the 
pressure term as 

p dx dx^p) ^ ^ 

In the finite diflference scheme of ZEUS-MP, this gives 

v^,i{t + At)-v^,i{t) ^ A f I I \ 

At Axi \pi{t) Pi-i{t)J ^ ' 

Applying this method of calculating pressure forces to 
the initial SLS setup gives the analytically expected result 
in all but the central 2"^ cells. 

We also modified the ZEUS-MP gravity solver to im- 
prove its accuracy. We added the octupole term to the 
multipole expansion used in the calculation of the grav- 
itational potential boundary values and we modified the 
parallelizaton scheme for this routine. While the inclu- 
sion of an octupole term in the boundary value routines is 
not strictly necessary for a spherically symmetric collapse 
calculation, we included it as a precautionary measure in- 
tended to more accurately capture any unexpected devi- 
ations from spherical symmetry. With these changes, the 
gravity solver was found to reproduce the expected poten- 
tial of the singular logatrope with an error of only 0.1% in 
the innermost 2^ cells, falling rapidly to 3 x 10""*% at the 
edge of the sphere. 

3.4. Setting p{r) on the Grid 
3.4.1. Density Within the Core 
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The ZEUS-MP gravity solver demands that, when plac- 
ing a splierical distribution of mass on a Cartesian grid, 
scrupulous attention be paid to the distribution of mass 
between the sphere and the edge of the computational 
grid. The corners of the grid are especially important 
because any significant mass located there will strongly 
distort the spherical symmetry of the potential. To avoid 
this problem, we surrounded our spheres with a medium 
whose density was 1000 times less than that at the sphere's 
edge. This method produced highly spherically symmetric 
potentials. While this density contrast is greater than is 
expected near the edges of real molecular cloud cores, tests 
with lower values showed that all density contrasts > 30 
gave similar results. 

Because the multigrid gravity solver fails to produce ac- 
curate results in the region of a sharp density disconti- 
nuity, the jump in density at the edge of the sphere was 
smoothed over several cells using a Gaussian function. We 
also found that it was helpful to situate the edge of the 
sphere slightly away from the grid boundary. For a grid 
with dimensions Ixlxl, the initial diameter of the sphere 
was no more than 0.91. 

In order to lay down a singular density profile on a 
Cartesian grid, care must be taken in handling the cell 
containing the singularity. Our approach to was to numer- 
ically integrate equation (7) to determine the mass that 
belonged in the central, singularity-containing cell, and to 
divide by the cell's volume to determine its density. 

3.4.2. Density in the External Medium 

Initially, we elected to treat the low-density outer 
medium in a manner similar to that used by FC93. In 
that method, the outer medium was made isothermal and 
its sound speed was set to the value required to produce 
the desired truncation pressure at the edge of the sphere. 
Experiments showed that, while this medium remained at 
a constant pressure for long periods of time during the 
collapse, it would eventually develop unacceptably large 
fluctuations. 

A successful treatment of the outer medium was devel- 
oped by following the method of Boss and Black (1982), 
wherein a constant pressure the pressure required to 
truncate the initial c;ore is maintained on the surface of 
the core throughout collapse. The boundary condition 
is strictly outflow, allowing material to leave the core as 
needed but otherwise conserving the mass initially placed 
within the boundary. The velocity of material outside the 
core is zeroed. The effect is to simulate the collapse of a 
core which has only a finite reservoir of mass from which to 
draw. This is a likely scenario for the formation of a star 
in an isolated core, and is also consistent with observations 
of clustered star formation in p Oph, which indicate that 
each accreting object can access only a finite amount of 
material (Motte, Andre, & Neri 1998). 

3.5. A Note on Time Scales 

In accordance with the convention of Foster & Cheva- 
lier (1993) we define t = to be the moment during the 
collapse of the nonsingular logatropic sphere at which the 
density profile becomes singular. Thus, adjustments from 
the nonsingular to the singular configurations occur at 
times t <0. The same convention is applied to the singu- 
lar logatrope: t = marks the beginning of the singular 



collapse. However, in order to facilitate the plotting of 
some of the results of the nonsingular collapse on log-log 
scales, we define an alternative time scale, T = t+tsing > 
where tsing is the elapsed time between the beginning of 
the simulation and the formation of the singularity. Thus, 
T = marks the beginning of the simulation for both the 
singular and nonsingular cores. Consequently T = t for 
the singular core, since it begins from a singular density 
profile (i.e. tging = for the singular core). 

For the purposes of comparison with observations, we 
have chosen to scale most of our results using fiducial val- 
ues of the truncation pressure, Pg, and central temper- 
ature, Tc- In order to facilitate reinterpretation of our 
results for other values of Pg and T^, wc have also pro- 
vided dimensionless versions of the key results. In their 
dimensionless forms, densities are expressed in terms of 
the initial central density, Pc, speeds in terms of the ini- 
tial central sound speed, ac, radii in terms of ro, masses 
in terms of the total mass of the core, Mtot, and times 
in terms of the mean free-fall time of the core, tg. All of 
our results can be rescaled to other values of Pg and Tc 
by recalculating each of these quantities using the data in 
Table 1 and the equations in the Appendix of MP96. 



3.6. The Sink Cell 

The presence of a singularity in the density profile of ma- 
terial on the grid can lead to two primary problems. The 
first is that infalling material may shock as it encounters 
the singularity (i.e. the protostellar object). The second 
is that the ever-increasing densities and velocities in the 
region around the singularity inevitably cause a crippling 
trend toward tiny numerical time steps. In order to avoid 
these problems, we employed a "sink cell" approach, mod- 
eled after that of Boss and Black (1982). The use of a 
sink cell effectively isolates the singularity and the details 
of the flow around it from the flow on the rest of the grid. 

In our method, the innermost cube of 3 x 3 x 3 cells 
centered on the singularity were blocked off as a collec- 
tive "sink cell". Within this region, two types of mass 
are counted. The first is mass which is still on the com- 
putational grid and which therefore undergoes the usual 
hydrodynamic and gravitational interactions. The second 
type of mass within the sink is that which has passed into 
a central point mass; this mass does not interact hydrody- 
namically with the rest of the grid, but continues to exert 
its gravitational infiuence via a point mass potential. The 
density of grid material within the sink is set to a spatially 
uniform value, which is reset at each time step to the min- 
imum of the densities of the face-sharing neighbour cells 
(the results are insensitive to the exact density in this re- 
gion). We assume that the mass within the sink can be 
considered to have condensed onto the central protostellar 
object (Boss and Black 1982). Thus, any mass falling into 
the sink in excess of the spatially uniform density therein 
is subtracted from the grid and added to the central point 
mass at the sink's centre. Figure 2 shows a schematic di- 
agram of our sink cell. 

Because the density of grid material within the sink (i.e. 
excluding the point mass) is not representative of the true 
density in that region, we determined the pressure gradi- 
ents between the cells in the sink and those immediately 
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Table 1 
Scaling Parameters 







1.34 


0.956 


2.21 


0.959 


3.26 


0.956 



Dimensionless scaling parameters for the three nonsingular logatropes. These values allow for rescaling to values of 
and Tc other than those used herein (see text). 
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Fig. 2. — 2-D slice through the midplanc of the sink cell. The shaded cells are inside the sink, which is a 3 X 3 X 3 cube of cells. All of the 
shaded cells have the same density, which is equal to the minimum of the densities of their face-sharing neighbours. The location of the point 
mass is indicated. The thick line indicates the border of the sink, where both the gravitational potential gradients and the pressure gradients 
are determined by linear extrapolation from the adjacent cells. 
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outside by linear extrapolation using the densities in the 
adjacent two cells. The gravitational potential gradients 
(due to grid matter) around the edge of the sink were de- 
termined in the same way. This was necessary because 
the multigrid solution for the potential falters at the edge 
of the sink, where the density makes the abrupt change 
from a singular to a flat-topped proflle. At late times, the 
dominant gravitational link between the sink and the rest 
of the grid is through the point-mass potential, so the ex- 
trapolation of the potential gradients of material on the 
grid represents a relatively small correction. 

The sink cell can be activated at any time without sig- 
nificant effect on the collapse. Figure 3 illustrates this 
point: it shows the mass accretion profiles for three sim- 
ulations of the collapse of a 1 Mq nonsingular logatropic 
sphere, differing only in the presence or absence of a sink 
cell and its time of activation. The difference in the mass 
contained within the central 3"'' cells between simulations 
with and without the sink cell never exceeds 5%, and is 
< 0.1% for t > O.GGfff. The effect on the velocities near 
the centre is equivalently small. Both the density and the 
velocity far from the centre are essentially unaffected by 
the use of a sink. The small difference visible in Figure 
3 between the simulations in which the sink is activated 
before the singularity forms (T = 0.26ffi) and just after (T 
= 1.14taf) is attributed to the development of supersonic 
inflow in the region of the sink just before the formation 
of the density singularity. The supersonic inflow isolates 
the sink hydrodynamically from the flow on the rest of the 
grid, essentially reducing its interaction with the rest of 
the grid to a gravitational one. Thus, if the sink is turned 
on before this supersonic inflow develops, its effect on the 
flow near the centre of the core is felt more strongly. 

The use of the sink cell alleviates the problem of shrink- 
ing time steps, extending the time to which the simulation 
can be followed. Additional gains in speed can be made 
by using lower resolution simulations. We use the higher 
resolution of 255^ cells in two situations: for t < in the 
nonsingular collapse, to monitor the development of the 
singularity; and for times prior to the departure of the ex- 
pansion wave from the core in the singular collapse, for 
the purpose of comparison with the analytics. After these 
times, however, our primary goal in both types of simula- 
tion is to track the mass accretion history of the central 
object, so we switch to the lower resolution of 127^. 

After the singular density profile has been achieved, 
there is little difference between the accretion profiles for 
simulations conducted at 127^ and 255^ cells (see Fig. 4). 
When the simulation begins (T = in Fig. 3), there is a 
large difference in the mass contained within the two sinks 
simply because the sink in the higher resolution encom- 
passes a smaller physical volume, and therefore contains 
less mass. After 2.0ts, however, there is only a 3 x 10~^% 
difference between the sink masses at the two resolutions. 
This is because, by this late time, the accreted mass is de- 
termined almost entirely by the mass of the point mass 
(whose dependence on resolution decreases with time). 
Hence, to follow Macc{t) to later times after singularity 
formation, we arc justified in switching to the lower reso- 
lution to benefit from the improved speed. 

3.7. The Jeans Criterion 



Truelove et al. (1997) established a criterion for ensuring 
that artificial fragmentation does not arise due to Jeans in- 
stability in hydrodynamic simulations. They determined 
that artificial fragmentation will not occur if the ratio of 
the grid spacing to the Jeans length, which they call the 
Jeans number, J = Ax/Xj, is kept below 0.25. We adopt 
this same implementation of the numerical Jeans criterion, 
with one difference: for the logatrope, the Jeans length be- 
comes 

/ TT dP UAPc 

Using this version of Aj, we check the entire grid every few 
time steps to ensure that J < 0.25. 

In all of our simulations, the presence of a central den- 
sity singularity inevitably gives rise to violations of the 
numerical Jeans condition. These violations occur only in 
the immediate vicinity of the singularity, where the inflow 
is supersonic. This supersonic inflow effectively isolates 
the Jeans-violating cells from the flow on the rest of the 
grid. For this reason, and because we observe no frag- 
mentation or significant asymmetry anywhere on the grid, 
we conclude that these violations of the numerical Jeans 
condition near the central cell can be disregarded. 

3.8. Testing 

We tested the code in three ways: we compared test 
problems to known results, we examined the quality of 
equilibria generated, and we checked that our simulation 
procedure could reproduce the theoretical solution for a 
singular logatrope. The first two of these tests are dis- 
cussed here, and the results of the third can be found in 
§4. 

In order to establish that the code could reproduce 
known results when configured for our systems, we ran 
several of the problems in the NCSA test suite. We were 
able to exactly reproduce published results for the 1-D 
Sod shock tube (Stone & Norman 1992a). Further, our 
simulations of the gravitational collapse of an initially hy- 
drostatic, uniform-density, pressure-free sphere on a 3-D 
Cartesian grid showed excellent agreement with the ana- 
lytic results of Hunter (1962). 

Tests showed that the code and our routines could estab- 
lish and maintain stable logatropic equilibria. Recall that 
we expect nonsingular cores of radius R/ro < 24.37 to ex- 
hibit stable equilibria (eq. [8]). For a significantly subcrit- 
ical core {R/tq = 1.34 or M = IMq, for our usual choice 
of central temperature and truncation pressure) we found 
that equilibrium could be maintained for many free-fall 
times. After being initialized in hydrostatic equilibrium, 
this core oscillated about that equilibrium with velocities 
never surpassing 3 x 10^^ times the local sound speed for 
periods of 20 free-fall times (after which we halted the 
calculation). Because we are restricted to using uniform 
grids, the larger the radius, R/r^, of the simulated sphere, 
the poorer is our ability to resolve ro, and hence to dis- 
tinguish numerically between singular (unstable) and non- 
singular (stable) cores. For this reason, we are unable to 
obtain stable equilibria for cores with R/ro > 4.5 (i.e. 
> 9Mq for our standard truncation pressure and central 
temperature), as these cores behave essentially indistin- 
guishably from singular cores (even at the higher resolu- 
tion of 255^ cells). Our study of nonsingular cores is there- 
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Fig. 3. — Accreted central mass, Mace, in the sink as a function of time for a 1 Mq nonsingular collapse. Mass is expressed as a fra<;tion 
of the total mass of the core, and time as a fraction of the free-fall time. The solid line represents a simulation with no sink cell. The others 
represent simulations with the sink activated at T = 0.26fff (dot-dashed line) and T = l.Utg (dashed line). The behaviour of the three 
solutions around the time of singularity formation (T = 1.09ffj = tging) is shown in the inset. 




Fig. 4. — Accreted mass within the sink, Mace as a function of time for 1 Mq nonsingular collapses at high (255"^ cells, dashed line) and 
low (127^ cells, solid line) resolutions. The vertical dotted line indicates the time at which the singular density profile is achieved. 
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fore confined to those with truncation radii R/rQ < 3.26 
(corresponding to masses < 5Mq), for which good initial 
equilibria can be obtained and vq can be properly resolved 
over 30 cells or more. To initiate collapse of the nonsin- 
gular logatropc, wc increased the density throughout the 
core by 5% above equilibrium values. 

Establishing and maintaining the unstable equilibrium 
of the singular logatrope is considerably more difficult. As 
expected, almost any small perturbation, including those 
inherent in the integration algorithms, is sufficient to cause 
it to collapse. Our standard procedure for the singular col- 
lapses was to set them up in the best numerical equilibrium 
we could achieve, and then initiate the collapse in a rela- 
tively controlled manner by increasing the density of the 
central cell by a few percent (the results are insensitive to 
the exact value of the density enhancement). 

3.9. Physical Core Conditions 

Given a value for A, the mass of a logatropic sphere is 
wholly determined by its radius, R, central temperature, 
Tc, and truncation pressure, Ps- In all of our simulations, 
we used A = 0.2 in the logatropic EOS. We take as our 
"standard" central temperature and truncation pressure 
the values (T^, P,) = (10 K, 1.3 x lO^ks ctr-^ K), chosen 
to be consistent with those used by Shu (1977) and MP97. 
These physical conditions are typical of those found in re- 
gions of relatively isolated star formation. For example, by 
fitting a Bonnor-Ebert sphere to the density structure of 
the isolated dark cloud Barnard 68, Alves, Lada, & Lada 
(2001) derived a surface pressure of 1.8 x 10^A;b cm~^ K 
and a temperature of 16 K. 

We emphasize that the Tc and Pg used herein are not 
representative of conditions in regions of clustered star for- 
mation, where the surface pressures can be one or two or- 
ders of magnitude larger (see §7). In order to facilitate the 
extrapolation of our results to different Ps and Tc values, 
we have included dimensionless versions of the primary 
results. 

Wc have simulated the collapse of nonsingular logat- 
ropes truncated at radii R/tq = 1.34, 2.21, and 3.26. Sub- 
ject to our standard Tc and Pg, these models have masses 
of 1, 2.5, and 5 M0 and scale radii of ro = 0.071,0.065, 
and 0.060 pc, respectively. The higher densities in more 
massive cores lead to greater optical depths, and hence the 
need to account for radiative effects. However, as our goal 
is to study the purely dynamical aspects of gravitational 
collapse when radiation pressure is not yet dominant, our 
models do not include radiative effects. Future simulations 
of the collapse of higher mass cores will have to account 
for such radiative effects. 

The singular logatrope considered herein is the pure 
p oc r^^ case of equation (7), corresponding to C = 
and Too = 6.67 x 10^'', in the terminology of MP97. By 
imposing a truncation radius and our standard physical 
conditions on the singular logatrope, we scale it such that 
its total mass is 1 Mq. 

3.10. From 3-D to 1-D 

In this paper, we often represent physical quantities 
(density and infall speed) as a function of radius within 
a core. To convert between a given quantity represented 
on a 3-D Cartesian grid and the same quantity expressed 



as a function of radius alone, we simply binned the data 
on the Cartesian grid into spherical shells, each one pixel 
thick and centred on the core's centre. We then averaged 
the given quantity over each shell to plot the quantity as 
a function of radius alone. 

4. SINGULAR COLLAPSE 

In this section we describe the results of simulations of 
the collapse of singular logatropes. Our objectives in sim- 
ulating the singular collapse are (a) to ensure that our 
method for simulating a collapse in 3-D accurately repro- 
duces the self-similar solution of MP97 and (b) to investi- 
gate the accretion history of a singular collapse in which 
the available mass reservoir is finite. 

Figure 5 shows the time evolution of the radial infall 
speed for a singular collapse (scaled such that the total 
mass of the core is 1 Mq). The solid lines in the figure are 
the analytic solutions of MP97 for the SLS collapse. Each 
one represents the infall speed of the gas as a function of 
radius at a given instant. 

In general, the core collapses in the expected inside-out 
manner described by MP97. Due to the initial 25% density 
enhancement applied to the central cell, the core is initially 
slightly out of hydrostatic balance and begins to collapse 
slowly at all radii. This effect can be seen in the figure 
as the low-speed deviations of the simulation data from 
the analytic curves, just past the expansion wave front at 
each time interval shown. At early times the expansion 
wave, which should be spherical, is poorly resolved on the 
rectangular grid, so we do not expect good agreement be- 
tween simulation and theory. By 0.5 fff, the expansion 
wave is easily visible over the small inward motion of the 
whole core. By 1.0 iff, the simulation shows good agree- 
ment with the theory. The density of the collapsing core, 
shown in Figure 6, also shows good agreement with the 
analytics of MP97. 

Once the expansion wave leaves the core, wc can no 
longer follow its behaviour, but we can continue to track 
the motions of material inside the core. As the 2.5 is line 
in Figure 5 indicates, the collapsing material within the 
truncation radius continues to closely follow the analytic 
predictions of MP97. 

Equation (10) tells us that the mass of the central ob- 
ject at r = in a collapsing SLS should increase as t^. 
We can measiirc an approximately similar quantity in our 
simulations by recording the accretion history, Macc(i)) 
of the sink cell. Figure 7 shows a plot of the accretion 
history of the 1 singular core. Also shown is a three- 
parameter fit to Mace = Co (i + Ci )^ + ^"2 , where the fitted 
parameters arc the Ci. Here, Co is a scaling factor play- 
ing the same role as the constant (APc47rG)'^/^ in equa- 
tion (10). The constant Ci accounts for the fact that, at 
t = Q, the simulation is already slightly out of equilib- 
rium, due to the applied density enhancement. Thus, the 
equilibrium state is effectively pushed to a time — Ci. For 
a wide variety of parameters (resolutions, sink cell acti- 
vation times, initial density enhancements, etc.) we find 
that Ci ~ tff/5. The last constant, C2, allows for the 
fact that some mass is already present within the sink 
when it is activated. The best-fit values of the Cj are 
(Co, Ci, C2) = (2.75 X 10-3, 0.210, 1.03 x lO'^). 

As predicted by MP97, the accreted mass in the SLS 
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Fig. 6. — Density, p, as a function of radius, r, within a collapsing singular core for times prior to l.Stff . The core is scaled such that 
-Mtot = 1 Mq. The meanings of symbols and lines are as in Fig. 5. Only that portion of the core in which the p(r) profiles are distinguishable 
from each other is shown. 
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collapse increases as prior to 1.8%, at which time the 
expansion wave leaves the core. After this time, the self- 
similarity of the solution is broken and we expect the sim- 
ulation to diverge from the self-similar SLS collapse so- 
lution. Because our bomidary condition is such that the 
accreting central object has a finite mass reservoir from 
which to draw (see §3.4), we can predict an eventual slow- 
down in the accretion rate. Once the expansion wave 
leaves the core, it is no longer setting new material into 
collapse, so there is no new flux of inward-moving material 
to maintain the flow. 

As shown in Figure 8, accretion onto the central object 
continues as Mace until t ~ 3 -Its, at which time 35% 
of the total initial mass of the core is bound up within the 
central point mass. Prom this point on, the accretion rate 
declines slowly until wc stop the simulation at 7.2ts, when 
90% of the mass of the core has been accreted. MP97 pre- 
dicted that 100% of the mass in an SLS collapse would 
be accreted by 4.3%, but this prediction is based on the 
assumption that the self-similar expansion wave solution 
holds over the entire collapse. Wc have simulated a core 
with a finite size, so the simulation ceases to be self-similar 
once the expansion wave has left the core. As our simula- 
tions show, the accretion timescalc must then increase, as 
no new material is being set to collapse by the passage of 
the expansion wave. 

5. NONSINGULAR COLLAPSE 

The evolution of the 1 Mq nonsingular core prior to the 

formation of the singular density profile is documented in 
Figures 9 and 10. The nonsingular collapse proceeds in a 
manner qualitatively similar to the Bonnor-Ebert collapse 
of FC93: collapse begins at all radii simultaneously (due to 
the initial applied overdensity), but the peak in the veloc- 



ity profile migrates inward slowly. The collapse sequence 
is somewhat different from that seen in the Bonnor-Ebert 
sphere. FC93 found that the inner, initially flat-topped 
region of the Bonnor-Ebert sphere collapses to an r^^ den- 
sity profile, thereby matching the power-law density profile 
of the rest of the sphere. We find that the initially flat- 
topped inner region of the nonsingular logatropic sphere 
collapses directly to an approximately r^^^^ density pro- 
file. Thus, at the moment when the inner core becomes 
truly singular — the moment we label as t = — the whole 
core has a dual power-law density structure: r"^^"^ for 
r < To and the initial r^^ for r > ro. Recall from §2.3 
that an r~^/^ profile is the form of the density expected 
at small radii in the collapse of the singular logatropic 
sphere. This would seem to indicate that we are seeing a 
fairly rapid transition to an SLS-like collapse. Note, how- 
ever, that it is only after this singular density profile is 
fully established in the inner core that we see a marked in- 
crease in the central accretion rate (i.e. the transition from 
the collapse phase to the accretion phase) . This transition 
will be discussed further in §6. 

To facilitate comparison between them. Figure 7 shows 
the mass accretion profile for the 1 Mq nonsingular col- 
lapse in addition to that of the singular collapse. Prior 
to the formation of the singularity, the rate of accretion 
onto the central mass increases slowly. Consequently, very 
little of the mass is accreted during this time: by t = 0, 
the sink contains only 0.4% of the total mass of the core. 
Just prior to the emergence of the singular density, the 
accretion rate increases markedly. It peaks just after the 
singularity forms, then enters a period of relatively steady 
accretion until about half of the total mass has been ac- 
creted. Thereafter, the accretion rate begins a slow decline 
which continues for the remainder of the collapse (see §6 
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Fig. 9. — Density evolution of the 1 Mq (R/ro = 1.34) nonsingular core. The singular density profile is reached at t = 0. The thick straight 
line indicates p oc and the thinner straight line indicates a p oc r~^/^ profile. The solid curve (t = — 1.09fff ) represents the initial density 
profile. 
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Table 2 
Accretion Histories 



t/ts 


t 


Macc/Mtot 


dM^c/dt 




(10^ yr) 




(10-6 Mo yr-i) 


1M0 


-0.504 


-0.244 


2.36 X 10-4 


7.22 X 10-* 


0.00 


0.0 


4.7 X 10-3 


0.167 


0.500 


0.242 


0.164 


0.702 


1.00 


0.484 


0.340 


0.760 


1.50 


0.725 


0.522 


0.712 


2.00 


0.967 


0.670 


0.499 


2.50 


1.21 


0.768 


0.325 


3.00 


1.45 


0.831 


0.211 


3.50 


1.69 


0.873 


0.139 


4.00 


1.93 


0.901 


0.0964 


4.50 


2.18 


0.921 


0.0702 


5.00 


2.42 


0.936 


0.0522 


5.43 


2.62 


0.945 


0.0425 


2.5MQ 


-0.505 


-0.286 


2.78 xlO-* 


1.41 xlO-^ 


0.00 


0.00 


3.80 xlO-=^ 


0.238 


0.500 


0.284 


0.140 


1.47 


1.00 


0.567 


0.310 


1.56 


1.50 


0.850 


0.488 


1.50 


2.00 


1.13 


0.635 


1.07 


2.50 


1.42 


0.734 


0.698 


3.00 


1.70 


0.798 


0.451 


3.50 


1.98 


0.840 


0.294 


4.00 


2.27 


0.868 


0.204 


4.50 


2.55 


0.888 


0.148 


5.00 


2.83 


0.902 


0.112 


5.33 


3.02 


0.910 


0.0940 


5M0 


-0.500 


-0.319 


4.97 xlO-'^ 


3.29 xlO-* 


-0.253 


-0.161 


7.84 xlO-^ 


2.11 xlO-3 


0.00 


0.00 


1.00 xlO-3 


0.179 


0.250 


0.159 


0.0332 


1.70 


0.500 


0.318 


0.106 


2.51 


0.750 


0.478 


0.186 


2.58 


1.07 


0.685 


0.294 


2.72 



Accretion history of the 1, 2.5, and 5 Mq, nonsingular, A = 0.2 logatropic spheres with our standard physical conditions 
of Pg = 1.3 X lO^fcs K cm-3 and = 10 K. Recall that t = marks the time at which the singular density profile is 
established. 
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Fig. 10. — Evolution of the radial infall speed of the 1 Mq [R/ro = 1.34) nonsingular logatrope preceding the development of the singular 
density profile. 



for further discussion). Table 2 shows a representative set 
of Mace and M^cc values for the nonsingular collapse. 

6. HIGHER MASS COLLAPSES 

Figure 11 shows the time evolution of the accreted 
central mass, Macc(i)i and central mass accretion rate, 
Macc{t), for the collapse of nonsingular spheres with 
R/ro = 1.34, 2.21, 3.26. To situate these results in a phys- 
ical context. Figures 11a and lib are scaled to according 
to our standard Pg and Tc, which gives them masses of 1, 
2.5, and 5 Mq. Dimensionless versions of the same data 
are supplied in Figures 11c and lid to facilitate scaling to 
other combinations of Pg and Tc- The R/r^ = 3.26 lines 
are truncated shortly after the singular density profile is 
achieved due to a limitation of the simulation technique. 
For this core, switching to the lower grid resolution to gain 
the necessary speed improvement would lead to insufficient 
resolution of the scale radius, vq. Thus, the simulation 
ends when the time step becomes untenably small at the 
higher resolution. 

According to the singular logatropic collapse theory of 
MP96, more massive cores are characterized by higher in- 
fall rates, but longer overall accretion lifetimes. We find 
that the nonsingular cores exhibit similar behaviour. Con- 
trast this with the SIS model, in which all cores have the 
same accretion rate, regardless of their total mass. There- 
fore, the accretion lifetimes of higher mass protostars are 
proportionately longer than in the logatropic scenario. As 
can be seen in Figure lib, the central accretion rate (in 
physical units) increases with increasing total core mass. 
Further, as shown in Figures 11a and 11c, the 1 M0 core 
is the first to reach Macc/Mtoi = 0.9; that is, the more 
massive cores have longer accretion timescales. The trend 
seen in Figure 11 then follows logically: in order for the 
core with the lowest mass to be the first to accrete 90% 



of its mass, it must accrete a larger fraction of its total 
mass, per unit time, than the higher mass cores. Thus, its 
accretion rate expressed in dimensionless units, must be 
higher than those of the more massive cores. 

The accretion histories of all three nonsingular logat- 
ropic models are summarized in Table 2. 

As in the collapse of the Bonnor-Ebert sphere, in the 
higher mass collapses, the central mass of the nonsingular 
logatropes accretes very slowly, prior to the formation of 
the singularity. In all three cases, the elapsed time be- 
tween the initiation of collapse and the formation of the 
density singularity is tging ~ 1.44t, where r = ro/cc is the 
sound crossing time of the region interior to tq. This is 
also just larger than the mean free-fall time of the region 
interior to tq. Recall that the initial density profile in the 
outer regions (r > ro) of the nonsingular logatrope closely 
follows the power law profile of the singular logatrope (see 
Fig. 1). Thus, it is reasonable that the time taken for 
the global readjustment to power law density behaviour 
should be of order the free-fall time of the initially fiat- 
topped region, r < vq. This information is summarized in 
Table 3. 

For all three masses studied, immediately following the 

formation of the singularity, the central mass accretion 
rate increases abruptly. The inflow then moves from the 
collapse phase into an SLS-like accretion phase similar to 
that seen in Figure 8. 

Simulations of the collapse of a critical Bonnor-Ebert 
sphere performed by FC93 showed that, at the moment 
the density profile becomes singular, 44% of the mass is 
in supersonic infall with Macli muubers as high as ~ 3.25. 
MP96 predicted that, because the logatropic EOS is softer 
than the isothermal EOS, a nonsingular logatropic collapse 
might not give rise to such high infall velocities. Our re- 
sults are consistent with this prediction. We have found 
that, while highly subcritical logatropes are characterized 
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Fig. 11. — Time evolution of Macc(t) and Macc(t) for the collapse of nonsingular logatropic spheres with R/ro = 1.34 (solid line), 2.21 
(dashed line), and 3.26 (dotted line). Panels (a) and (b) show the data scaled according to our standard paxameters such that the models 
have total masses of 1, 2.5, and 5 Mq. Panels (c) and (d) show the same data in dimensionless form. The vertical dashed lines in panel (b) 
indicate the suggested transitions between the PPC, Class 0, and Class I stages for the 1 model (see text). 



Table 3 

Times of Singularity Formation 



Mass 


r 


ising 


ising 


(Mo) 


(10^ yr) 


{r) 


(te('^o)) 


1 


4.84 


1.44 


1.29 


2.5 


5.67 


1.45 


1.35 


5.0 


6.36 


1.44 


1.40 



Elapsed time between initiation of collapse and singularity formation, fging) for the three nonsingular logatropic spheres 
subject to our standard physical conditions. 
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by supersonic infall at the moment of singularity forma- 
tion, the trend is toward a gentler collapse as we approach 
criticality (i.e. as we increase the radius of the core). 

Table 4 lists the percentage of the mass of each simulated 
logatropic core which is inf ailing supcrsonically at t — 0. 
The clear trend is that, as the radius of the core increases, 
both the fractional amount of mass which is in supersonic 
infall and the Mach number of that infall decrease. Ex- 
trapolating the trend, it seems likely that little or none of 
the mass of a critical R/vq = 24.37 logatropc would be in 
supersonic infall at f = 0, in accord with the prediction of 
MP96. Simulations of the collapse of a critical core will be 
needed to confirm this prediction. 

7. DISCUSSION 

As star-forming cores evolve from Class to Class IV, 
their protostellar mass accretion rate may vary by ~ 2 
orders of magnitude (see Andre, Ward- Thompson, & Bar- 
sony 2000, and references therein). It should be noted 
that estimates of the accretion rate in young stellar ob- 
jects are not made by direct measure, but rather by in- 
ference from a combination of models and other data, as 
will be discussed shortly. In contrast to the observations, 
the standard SIS model predicts a time-invariant proto- 
stellar mass accretion rate. Foster & Chevalier (1993) and 
Boss and Black (1982) have obtained variable mass ac- 
cretion rates for isothermal collapse models which draw 
mass from finite reservoirs, but these models show a ten- 
dency to evolve toward the constant mass accretion rate 
predicted by the SIS model. Henriksen, Andre, & Bon- 
temps (1997) and Caselli and Myers (1995) have proposed 
alternative models for protostellar collapse which predict 
variable accretion rates. Like the logatrope, these mod- 
els are phcnomenologically based, designed to match the 
observed properties of star-forming cores (density profiles, 
nonthermal linewidths, etc.) These models, however, have 
many more free parameters than the logatrope, which has 
only one: the constant A in equation (3). As shown in the 
previous section, the logatrope also predicts a highly vari- 
able accretion rate. We must ask, then, how closely the 
logatropic accretion profile corresponds to available obser- 
vational data. 

Bontemps et al. (1996) surveyed 36 Class I and 9 Class 
protostars in the p Oph, Taurus- Auriga, and Perseus star- 
forming regions and found that the outflow momentum 
flux, Fco declines from ~ 10~^ M0 km s~^yr~^ in the 
Class sources to ^ 2 x 10~^ Mq km s^^yr^-'^ in the Class 
I sources. Their suggestion was that this decrease in Fqo 
and the accompanying decrease in the mass ejection rate 
of the protostellar wind, M^^, were attributable to a de- 
crease in the protostellar mass accretion rate. Mace, from 
~ IQ-^ Mq yr-i to ~ 2 x 10"'^ Mg yr'^. Clearly the 
largest of these accretion rates is significantly higher than 
those shown in Figure 11, but we note that there is sig- 
nificant room for movement in both the observations and 
the scaling of our simulations. 

We have scaled our accretion profiles to central temper- 
atures and truncation pressures which are characteristic 
of regions of isolated star formation. If we instead choose 
Tc and Pg in keeping with observations of regions such as 
p Oph, we can make up the difference between our Mgcc 
values and those suggested by Andre, Ward- Thompson, & 



Barsony (2000). Johnstone et al. (2000) fitted 55 molec- 
ular cloud cores in p Oph to Bonnor-Ebert spheres and 
thereby made estimates of the surface pressure on each 
core. The typical truncation pressures found in that study 
lay in the range Pg = 10® — lO^fc^ cm^^ K. These arc be- 
tween 1 and 2 orders of magnitude greater than the fiducial 
truncation pressure used herein (Pg = 1.3 x lO^ks cm^ K). 
The relationship between the total mass of a nonsingular 
logatropic core and its truncation pressure takes the form 

— 1/2 

Mtot Ps ' ■ Hence, increasing the truncation pressure 
by a factor of 100 would reduce the mass of a given core by 
a factor of 10. This would also bring the mass of the crit- 
ical A = 0.2 logatrope down from 92 Mq — typical of an 
entire molecular cloud clump — to 9.2 Mq, which is typ- 
ical of a molecular cloud core situated within a clump. 
Hence, our models would then represent cores of compara- 
tively low mass. Wc would then have to look at logatropic 
spheres closer to the critical radius in order to find cores 
of mass ~ 1 Mq, in which case the accretion rates ought 
to be considerably higher than those listed in Table 2. 

The method used by Johnstone et al. (2000) to extract 
truncation pressures was, however, dependent on the as- 
sumption that all of the cores could be well-fit by Bonnor- 
Ebert spheres, and hence the truncation pressures ex- 
tracted cannot necessarily be considered reflective of those 
that would be required to truncate logatropic spheres. The 
dimensionless results of Figures 11a and lib are meant to 
be rescalable, should future EOS-independent measures of 
Pg become available. 

Disregarding the absolute numbers for a moment (which 
are, as we have emphasized, subject to considerable un- 
certainty), we note that the trend in our logatropic Mace 
bears significant resemblances to the data. Estimates place 
the transition from the Class to the Class I protostellar 
phase at the point when the accreted mass is approxi- 
mately equal to that remaining in the protostellar enve- 
lope. Mace — Menv (Andre, Ward-Thompson, & Barsony 
2000). As shown in Table 3, the nonsingular logatropic 
collapse behaves in just this way once about half of the 
mass of the core has been accreted, the accretion rate, 
Mace, begins to fall off. This behaviour arises because the 
accreting central object has only a finite mass reservoir 
from which to draw (see discussion in § 4). We suggest 
that one interpretation of the mass accretion profiles in 
Figure 11 is that they describe the transition from Class 
to Class I objects. If we take the 1 Mq case as an exam- 
ple, we sec that, at times t < 0, there is little accretion 
activity, so this may correspond to a very early Class or 
even a preprotostellar stage. From t ~ to ~ 0.7 x 10® 
yr, at which time the mass of the protostar and the re- 
maining envelope are equal, the object undergoes a period 
of vigorous accretion in an SLS-like manner. During this 
interval, the accretion rate, while variable, stays relatively 
constant, never deviating by more than 6% from the aver- 
age of 7.4 x 10"'' Mq yr-^ 

If the relationship between the rate at which material is 
accreted onto a young protostar and the rate at which mass 
is ejected through bipolar jets or outfiows were known, it 
would be possible to extrapolate the expected outflow rate 
from this predicted accretion rate. Unfortunately, no di- 
rect measure of the infall/outflow relationship exists, so 
we must resort to making an educated guess based on the 
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Table 4 
Supersonic Infall Properties 



i?/ro 


% of Mass Infalling 


Highest f/i^sound 




Supersonically at t = 


att = 


1.34 


5.5 


5.7 


2.21 


2.7 


4.9 


3.26 


1.3 


3.2 



Percentage of core's total mass which is infalling supersonically at f = 0, when the singularity forms. Note the decreasing 
trend in both percentage of mass in supersonic infall and the maximum infall Mach number with increasing R/vq- 



available models. Following Andre, Ward- Thompson, & 
Barsony (2000) and Pelletier & Pudritz (1992), we assume 
that realistic jet models give Afjet/Macc ^ 0.1 — 0.3. With 
this assumption, we extrapolate that accretion at the rate 
of 7.4 X 10"'' Mq yr~^ would be sufficient to power an 
outflow with an average Mjet ~ 0.74 - 2.2 x 10"'' M© 
yr""'^. 

In Figure 11, we have indicated our proposed evolution- 
ary sequence for the 1 Mq nonsingular logatropic core from 
the preprotostellar to the Class I stage. The preproto- 
stellar core (PPC) stage begins with the core in its near- 
equilibrium hydrostatic state. The PPC stage ends when 
the flat-topped central region of the core has collapsed to a 
singular r~^^^ density profile, just at the onset of vigorous 
accretion. This marks the end of the collapse phase and 
the transition to the accretion phase. The Class stage 
is characterized by vigorous accretion throughout. In ac- 
cordance with the description of Andre, Ward-Thompson, 
& Barsony (2000), we mark the transition from the Class 
to the Class I stage when Mace = M^nv Once in the 
proposed Class I stage, the accretion rates of our simu- 
lated cores decline rapidly. This decline in the accretion 
rate correlates well with the observed absence of strong 
outflows in Class I cores. 

One of the appeals of the logatropic model is its rela- 
tively small spread in formation times for stars of a wide 
range of masses. The accretion timescales for the forma- 
tion of 1-5 M0 protostars from nonsingular logatropes are 
on the order of a few million years, for the combinations of 
R/ro,Ps, and Tc used herein. These accretion timescales 
agree well with the observed ages of young star clusters, 
such as the Orion Nebular Cluster (ONC), wherein the 
bulk of the stars formed within ~ 10^ yr (Hillenbrand 
& Carpenter 2000). We anticipate that the accretion 
timescale will be shorter in regions of higher pressure, as 
in regions of clustered star formation. 

More sophisticated models of collapsing molecular cloud 
cores will include the effects of rotation and magnetic 
fields. The former ensures that most of the material that 
accretes onto the protostar does so through a circumstel- 
lar accretion disk. The latter, by contrast, is expected 
to cushion the collapse and perhaps to reduce the infall 
speeds. Magnetic fields will also enhance the function of 
large "pseudo-disks" (Galli & Shu 1993a) through mag- 
netic focusing of infalling material. It is important to em- 
phasize that the MP97 models include a measure of mag- 
netic support, in addition to the turbulent support. The 
logatrope incorporates a mean magnetic field as a virial pa- 
rameter, but actual time-dependent simulations employing 



a full 3-D magnetic field structure are necessary before one 
can gain a complete picture of logatropic collapse. 

8. SUMMARY 

We have conducted 3-D numerical hydrodynamic sim- 
ulations of the collapse of singular and nonsingular loga- 
tropic spheres, starting from initial states of hydrostatic 
equilibrium. The primary conclusions of this research are 

as follows: 

1. We have extended the solution for the SLS collapse 
beyond the strict limit of applicability of the 

analytic expansion wave solution {t < l.Stff). 
For an SLS collapse which draws from a finite 
mass reservoir, we find that the mass of the 
central accreting object increases initially as 
t"*, in accordance with the theory of MP97. 
This behaviour continues until t = S.lis, when 
Macc/Mtot = 0.35. 

2. We have developed Macc(i) and Macc(i) profiles for 
the collapse of nonsingular logatropic spheres with 
R/ro = 1.34,2.21, and 3.26. These values oi R/ro 
correspond to 1, 2.5, and 5 Mq cores, for our choice 
of fiducial truncation pressure, Pg = 1.3 x lO^ks 
cm~^ K, and central temperature, Tc = 10 K. 
We emphasize the accretion results because they 
pertain directly to observations of the luminosity 
of outfiows associated with protostellar cores, and 
therefore serve as a useful tool for distinguishing 
between molecular cloud core collapse models. 
The nonsingular logatropes are all characterized 
by initially slow accretion, then a period of intense 
accretion immediately following the development 
of a singular density profile. The intense accretion 
continues until Afacc — 0.5Aftoi! after which 

the accretion rate declines gradually until the 
entire mass of the core is accreted. For the three 
nonsingular collapses, we find that the elapsed time 
between the initiation of collapse and the moment 
of singularity formation (i.e. the formation of 
the hydrostatic protostellar object) obeys the 
approximate scaling tging — 1.44ro/crc ~ tg{ro)- 
The values of Pg and Tc used here are typical of 
isolated star formation; rescaling to other pressures 
and temperatures would allow the reinterpretation 
of these results in terms of clustered star formation. 

3. We note that the form of Macc(i) is qualitatively 
suggestive of a transition from the Class to 
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Class I protostellar phases. Improved knowledge 
of the precise physical conditions (Tg and P,) in 
the region of specific star-forming cores will be 
required to enable quantitative verification of this 
hypothesis. 

4. The trend in the accretion profiles is toward 
dominance of the envelope-type (i.e. p oc r~^) 
accretion phase with increasing R/ro (see Fig. 11). 
Simulations of cores with higher values oi R/tq 
will be needed to confirm this trend. 

5. Our results indicate that, in keeping with 
expectations based on the softness of the logatropic 
EOS (MP96), the adjustment of the nonsingular 
logatropic sphere to a singular density profile 

occurs more gently than in similar isothermal 
collapses (FC93). Wc showed that, although 
the collapse of subcritical nonsingular logatropes 
gives rise to supersonic infall at the moment of 



singularity formation, the fraction of the infalling 
mass which is in supersonic infall and the Mach 
number of the supersonic fiow both decrease 
as i?/ro increases towards the critical value. 
Highly-resolved simulations of the collapse of the 
critical nonsingular logatrope will be required to 
verify that its approach to the singular density 
profile is entirely subsonic. 
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